# PoP - Policing Socio-Geographic Boundaries and Inequality
# script for creating Figure 4
# Figure 4 - Influence of Logged Crime on Logged Police Stops (Standardized), Conditional on Racial Boundary Status.

# load packages
suppressPackageStartupMessages(
  
  {
    library(dplyr)
    library(tidyverse)
    library(ggplot2)
    library(haven)
    library(readxl)
    library(readr)
    library(areal)
    library(car)
    library(estimatr)
    library(magrittr)
    library(texreg)
    library(sandwich)
    library(jtools)
    library(prediction)
  }
)

## Austin

load("aus_stops_final.RData")

# now run analyses

# log and +1 to variables (pop already logged)
aus_stops_final = aus_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# now for new stop vars
aus_stops_final = aus_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
aus_stops_final$larrest_sd <- aus_stops_final$larrests - (mean(aus_stops_final$larrests)/sd(aus_stops_final$larrests))
aus_stops_final$lmisdemeanors_sd <- aus_stops_final$lmisdemeanors - (mean(aus_stops_final$lmisdemeanors)/sd(aus_stops_final$lmisdemeanors))
aus_stops_final$lfelonies_sd <- aus_stops_final$lfelonies - (mean(aus_stops_final$lfelonies)/sd(aus_stops_final$lfelonies))
aus_stops_final$lnonviolent_sd <- aus_stops_final$lnonviolent - (mean(aus_stops_final$lnonviolent)/sd(aus_stops_final$lnonviolent))
aus_stops_final$lviolent_sd <- aus_stops_final$lviolent - (mean(aus_stops_final$lviolent)/sd(aus_stops_final$lviolent))
aus_stops_final$lsociety_sd <- aus_stops_final$lsociety - (mean(aus_stops_final$lsociety)/sd(aus_stops_final$lsociety))
aus_stops_final$lperson_sd <- aus_stops_final$lperson - (mean(aus_stops_final$lperson)/sd(aus_stops_final$lperson))
aus_stops_final$lproperty_sd <- aus_stops_final$lproperty - (mean(aus_stops_final$lproperty)/sd(aus_stops_final$lproperty))

# create scaled DVs for stop vars
aus_stops_final$lall_stops_total_sd <- aus_stops_final$lall_stops_total - (mean(aus_stops_final$lall_stops_total)/sd(aus_stops_final$lall_stops_total))
aus_stops_final$lall_stops_black_sd <- aus_stops_final$lall_stops_black - (mean(aus_stops_final$lall_stops_black)/sd(aus_stops_final$lall_stops_black))
aus_stops_final$lall_stops_latino_sd <- aus_stops_final$lall_stops_latino - (mean(aus_stops_final$lall_stops_latino)/sd(aus_stops_final$lall_stops_latino))
aus_stops_final$lall_stops_white_sd <- aus_stops_final$lall_stops_white - (mean(aus_stops_final$lall_stops_white)/sd(aus_stops_final$lall_stops_white))
aus_stops_final$lall_stops_asian_sd <- aus_stops_final$lall_stops_asian - (mean(aus_stops_final$lall_stops_asian)/sd(aus_stops_final$lall_stops_asian))
aus_stops_final$lall_stops_nonwhite_sd <- aus_stops_final$lall_stops_nonwhite - (mean(aus_stops_final$lall_stops_nonwhite)/sd(aus_stops_final$lall_stops_nonwhite))


# create binned boundary measure
aus_stops_final <- aus_stops_final %>% mutate(boundary_quart = quantile(aus_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

# read in commercial density

# merge with commercial density data
load("aus_cd.RData")

head(aus_cd)

aus_cd_avg <- aus_cd %>%
  group_by(BLK_CODE) %>%
  summarize(total_jobs = mean(total_jobs, na.rm = TRUE))

aus_stops_final <- merge(aus_stops_final,aus_cd_avg, by = "BLK_CODE")

aus_stops_final$com_density <- aus_stops_final$total_jobs/aus_stops_final$pop

# merge with phys_bound
load("aus_fin_pb3.RData")

aus_fin_pb3 = aus_fin_pb3 %>% as.data.frame() %>% dplyr::select(-c(geometry))

aus_fin_pb3 = aus_fin_pb3 %>% dplyr::select(c(BLK_CODE, phys_bound))
aus_stops_final <- merge(aus_stops_final, aus_fin_pb3, by = "BLK_CODE")


###### race_blvXcrime plots #####
## Austin
mdata <- aus_stops_final %>%dplyr::select(lall_stops_total_sd,boundary_quart_dummy,
                                          lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                          hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                          com_density, phys_bound) %>% na.omit()


aus1 = lm_robust(lall_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + com_density + phys_bound, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,8.1,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  com_density = mean(com_density),
  phys_bound = mean(phys_bound)))

low_preds<- data.frame(cbind(prediction(aus1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,8.1,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  com_density = mean(com_density),
  phys_bound = mean(phys_bound)))

hi_preds<- data.frame(cbind(prediction(aus1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

aus_stops_preds <- preds %>% mutate(city = "Austin")


## Milwaukee
load("mil_stops_final.RData")


# log and +1 to variables (pop already logged)
mil_stops_final = mil_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# now for new stop vars
mil_stops_final = mil_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1),
         lped_stops_total = log(ped_stops_total + 1),
         lped_stops_black = log(ped_stops_black + 1),
         lped_stops_latino = log(ped_stops_latino + 1),
         lped_stops_white = log(ped_stops_white + 1),
         lped_stops_asian = log(ped_stops_asian + 1),
         lped_stops_nonwhite = log(ped_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
mil_stops_final$larrest_sd <- mil_stops_final$larrests - (mean(mil_stops_final$larrests)/sd(mil_stops_final$larrests))
mil_stops_final$lmisdemeanors_sd <- mil_stops_final$lmisdemeanors - (mean(mil_stops_final$lmisdemeanors)/sd(mil_stops_final$lmisdemeanors))
mil_stops_final$lfelonies_sd <- mil_stops_final$lfelonies - (mean(mil_stops_final$lfelonies)/sd(mil_stops_final$lfelonies))
mil_stops_final$lnonviolent_sd <- mil_stops_final$lnonviolent - (mean(mil_stops_final$lnonviolent)/sd(mil_stops_final$lnonviolent))
mil_stops_final$lviolent_sd <- mil_stops_final$lviolent - (mean(mil_stops_final$lviolent)/sd(mil_stops_final$lviolent))
mil_stops_final$lsociety_sd <- mil_stops_final$lsociety - (mean(mil_stops_final$lsociety)/sd(mil_stops_final$lsociety))
mil_stops_final$lperson_sd <- mil_stops_final$lperson - (mean(mil_stops_final$lperson)/sd(mil_stops_final$lperson))
mil_stops_final$lproperty_sd <- mil_stops_final$lproperty - (mean(mil_stops_final$lproperty)/sd(mil_stops_final$lproperty))

# create scaled DVs for stop vars
mil_stops_final$lall_stops_total_sd <- mil_stops_final$lall_stops_total - (mean(mil_stops_final$lall_stops_total)/sd(mil_stops_final$lall_stops_total))
mil_stops_final$lall_stops_black_sd <- mil_stops_final$lall_stops_black - (mean(mil_stops_final$lall_stops_black)/sd(mil_stops_final$lall_stops_black))
mil_stops_final$lall_stops_latino_sd <- mil_stops_final$lall_stops_latino - (mean(mil_stops_final$lall_stops_latino)/sd(mil_stops_final$lall_stops_latino))
mil_stops_final$lall_stops_white_sd <- mil_stops_final$lall_stops_white - (mean(mil_stops_final$lall_stops_white)/sd(mil_stops_final$lall_stops_white))
mil_stops_final$lall_stops_asian_sd <- mil_stops_final$lall_stops_asian - (mean(mil_stops_final$lall_stops_asian)/sd(mil_stops_final$lall_stops_asian))
mil_stops_final$lall_stops_nonwhite_sd <- mil_stops_final$lall_stops_nonwhite - (mean(mil_stops_final$lall_stops_nonwhite)/sd(mil_stops_final$lall_stops_nonwhite))


# now create binned racial blv measures
mil_stops_final <- mil_stops_final %>% mutate(boundary_quart = quantile(mil_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))
# merge with commercial density data
load("mil_cd.RData")

head(mil_cd)

mil_cd_avg <- mil_cd %>%
  group_by(BLK_CODE) %>%
  summarize(total_jobs = mean(total_jobs, na.rm = TRUE))

mil_stops_final <- merge(mil_stops_final,mil_cd_avg, by = "BLK_CODE")

mil_stops_final$com_density <- mil_stops_final$total_jobs/mil_stops_final$pop

# merge with phys_bound
load("mil_fin_pb3.RData")

mil_fin_pb3 = mil_fin_pb3 %>% as.data.frame() %>% dplyr::select(-c(geometry))

mil_fin_pb3 = mil_fin_pb3 %>% dplyr::select(c(BLK_CODE, phys_bound))
mil_stops_final <- merge(mil_stops_final, mil_fin_pb3, by = "BLK_CODE")


mdata <- mil_stops_final %>%dplyr::select(lall_stops_total_sd,boundary_quart_dummy,
                                          lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                          hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                          com_density, phys_bound) %>% na.omit()


mil1 = lm_robust(lall_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + com_density + phys_bound, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,7.5,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(mil1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,7.5,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(aus1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

mil_stops_preds <- preds %>% mutate(city = "Milwaukee")


## Chicago
load("chi_stops_final.RData")

# now create binned racial blv measures
chi_stops_final <- chi_stops_final %>% mutate(boundary_quart = quantile(chi_stops_final$p_race_white_blv, prob=.75, na.rm=TRUE),
                                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))

# now run analyses

# log and +1 to variables (pop already logged)
chi_stops_final = chi_stops_final %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(property_crime + 1),
         lviolentcrime = log(violent_crime + 1))

# now for new stop vars
chi_stops_final = chi_stops_final %>%
  mutate(lall_stops_total = log(all_stops_total + 1),
         lall_stops_black = log(all_stops_black + 1),
         lall_stops_latino = log(all_stops_latino + 1),
         lall_stops_white = log(all_stops_white + 1),
         lall_stops_asian = log(all_stops_asian + 1),
         lall_stops_nonwhite = log(all_stops_nonwhite + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
chi_stops_final$larrest_sd <- chi_stops_final$larrests - (mean(chi_stops_final$larrests)/sd(chi_stops_final$larrests))
chi_stops_final$lmisdemeanors_sd <- chi_stops_final$lmisdemeanors - (mean(chi_stops_final$lmisdemeanors)/sd(chi_stops_final$lmisdemeanors))
chi_stops_final$lfelonies_sd <- chi_stops_final$lfelonies - (mean(chi_stops_final$lfelonies)/sd(chi_stops_final$lfelonies))
chi_stops_final$lnonviolent_sd <- chi_stops_final$lnonviolent - (mean(chi_stops_final$lnonviolent)/sd(chi_stops_final$lnonviolent))
chi_stops_final$lviolent_sd <- chi_stops_final$lviolent - (mean(chi_stops_final$lviolent)/sd(chi_stops_final$lviolent))
chi_stops_final$lsociety_sd <- chi_stops_final$lsociety - (mean(chi_stops_final$lsociety)/sd(chi_stops_final$lsociety))
chi_stops_final$lperson_sd <- chi_stops_final$lperson - (mean(chi_stops_final$lperson)/sd(chi_stops_final$lperson))
chi_stops_final$lproperty_sd <- chi_stops_final$lproperty - (mean(chi_stops_final$lproperty)/sd(chi_stops_final$lproperty))

# create scaled DVs for stop vars
chi_stops_final$lall_stops_total_sd <- chi_stops_final$lall_stops_total - (mean(chi_stops_final$lall_stops_total)/sd(chi_stops_final$lall_stops_total))
chi_stops_final$lall_stops_black_sd <- chi_stops_final$lall_stops_black - (mean(chi_stops_final$lall_stops_black)/sd(chi_stops_final$lall_stops_black))
chi_stops_final$lall_stops_latino_sd <- chi_stops_final$lall_stops_latino - (mean(chi_stops_final$lall_stops_latino)/sd(chi_stops_final$lall_stops_latino))
chi_stops_final$lall_stops_white_sd <- chi_stops_final$lall_stops_white - (mean(chi_stops_final$lall_stops_white)/sd(chi_stops_final$lall_stops_white))
chi_stops_final$lall_stops_asian_sd <- chi_stops_final$lall_stops_asian - (mean(chi_stops_final$lall_stops_asian)/sd(chi_stops_final$lall_stops_asian))
chi_stops_final$lall_stops_nonwhite_sd <- chi_stops_final$lall_stops_nonwhite - (mean(chi_stops_final$lall_stops_nonwhite)/sd(chi_stops_final$lall_stops_nonwhite))


# merge with commercial density data
load("chi_cd.RData")

head(chi_cd)

chi_cd_avg <- chi_cd %>%
  group_by(BLK_CODE) %>%
  summarize(total_jobs = mean(total_jobs, na.rm = TRUE))

chi_stops_final <- merge(chi_stops_final,chi_cd_avg, by = "BLK_CODE")

chi_stops_final$com_density <- chi_stops_final$total_jobs/chi_stops_final$pop

# merge with phys_bound
load("chi_fin_pb3.RData")

chi_fin_pb3 = chi_fin_pb3 %>% as.data.frame() %>% dplyr::select(-c(geometry))

chi_fin_pb3 = chi_fin_pb3 %>% dplyr::select(c(BLK_CODE, phys_bound))
chi_stops_final <- merge(chi_stops_final, chi_fin_pb3, by = "BLK_CODE")


mdata <- chi_stops_final %>%dplyr::select(lall_stops_total_sd,boundary_quart_dummy,
                                          lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                          hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                          phys_bound, com_density) %>% na.omit()


chi1 = lm_robust(lall_stops_total_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + com_density + phys_bound, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,9.3,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(chi1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,9.3,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(aus1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

chi_stops_preds <- preds %>% mutate(city = "Chicago")


# graph out of sample predictions
aus_stops_preds <- aus_stops_preds %>% dplyr::select(lcrime,fitted.fit,
                                                     fitted.lwr,fitted.upr,boundary,city)

mil_stops_preds <- mil_stops_preds %>% dplyr::select(lcrime,fitted.fit,
                                                     fitted.lwr,fitted.upr,boundary,city)

chi_stops_preds <- chi_stops_preds %>% dplyr::select(lcrime,fitted.fit,
                                                     fitted.lwr,fitted.upr,boundary,city)



d <- rbind(aus_stops_preds, chi_stops_preds, mil_stops_preds)


# make plot
p <- ggplot(data=d, aes(x=lcrime, y=fitted.fit,group=boundary)) +
  facet_wrap(~ city, nrow = 3, scales = "free", strip.position = "right")+
  #geom_point(aes(shape = Boundary), color = "grey",alpha=.5)+
  geom_line(aes(y=fitted.fit, linetype=boundary, color = boundary))+
  geom_ribbon(aes(ymin=fitted.lwr, ymax=fitted.upr,linetype=boundary, color=boundary),alpha=0.1)+
  labs(y="Stops (Logged)",x="Crime (Logged)") +
  theme_bw(base_size=12,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA))+
  theme(strip.background = element_rect(fill="white"))

ggsave(plot = p, filename = "figure4_plots.png", 
       width = 8, height = 12, bg = 'white')


